#version 460

#define M_PI 3.1415926535897932384626433832795
#define R_MAX 0.01

layout(local_size_x = 1024, local_size_y = 1, local_size_z = 1) in;

uniform int tex_w, tex_h;
uniform ivec2 index_kernel_iv;

shared int WorkGroupSize;
shared int index_kernel_row, index_kernel_col;
shared int index_kernel;
shared vec3 pos_i_j;

layout(rgba32f, binding = 0) uniform image2D world_pos_map;
layout(rgba32f, binding = 1) uniform image2D kernel;

float fDiffuseProfile(float r, float A = 1.0, float s = 3.78);

void main() {
	int GlobalInvocationIndex = int(gl_WorkGroupID.x * gl_WorkGroupSize.x + gl_LocalInvocationID.x);
	if (GlobalInvocationIndex == 0)
	{
		WorkGroupSize = int(gl_WorkGroupSize.x * gl_WorkGroupSize.y * gl_WorkGroupSize.z);
		index_kernel_row = index_kernel_iv.y;
		index_kernel_col = index_kernel_iv.x;
		index_kernel = index_kernel_row * tex_w + index_kernel_col;
		pos_i_j = imageLoad(world_pos_map, ivec2(index_kernel_col, index_kernel_row)).xyz;
	}
	barrier();
	if (pos_i_j != vec3(0, 0, 0))
	{
		// Compute kernel at (i, j).
		for (int index_texel = GlobalInvocationIndex; index_texel < tex_h * tex_w; index_texel += WorkGroupSize)
		{
			int row = index_texel / tex_w;
			int col = index_texel % tex_w;
			vec3 pos_row_col = vec3(imageLoad(world_pos_map, ivec2(col, row)));
			float l = length(pos_i_j - pos_row_col);
			if (l <= R_MAX && l != 0)
			{
				float R = fDiffuseProfile(l, 1.0, 3.78);
				imageStore(kernel, ivec2(col, row), vec4(R, 0, 0, 0));
			}
			else
			{
				imageStore(kernel, ivec2(col, row), vec4(0, 0, 0, 0));
			}
		}
		barrier();
		for (int index_texel = GlobalInvocationIndex; index_texel < tex_h * tex_w; index_texel += WorkGroupSize)
		{
			int row = index_texel / tex_w;
			int col = index_texel % tex_w;
			imageStore(kernel, ivec2(col, row), vec4(imageLoad(kernel, ivec2(col, row)).x, 0, 0, 0));
		}
		barrier();
	}	
	else
	{
		if (GlobalInvocationIndex == 0)
		{
			imageStore(kernel, ivec2(0, 0), vec4(0, 0, 1, 0));
		}
		barrier();
	}
	barrier();
}

float fDiffuseProfile(float r, float A, float s) // A = 0.1, s = 3.78
{
	return A * s * ( ( exp(-s * r) + exp(-s * r / 3) ) / (8 * M_PI * r) );
}
